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CN ■ We study isotropic and anisotropic (Bianchi I) cosmologies in Palatini f{R) and 

^-,1 f{R, Rfj.uR^'^) theories of gravity and consider the existence of non-singular bouncing 
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Q^l gravitational solution to the big bang singularity and anisotropy problems of general 

\^ ' relativity without the need for exotic {w > \) sources of matter /energy. 
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I. INTRODUCTION 

Ever since its publication, Einstein's theory of general relativity (GR) has fascinated the- 
oretical physicists. Not only it is in excellent quantitative agreement with all observations 
[l[] (if a cosmological constant is included), but it also allows us to determine at which stage 
it should not be trusted. The existence of cosmological (big bang) and black hole singular- 
ities are clear symptoms that the theory is not complete. To overcome this drawback, it is 
common to argue that in such extreme scenarios quantum gravitational effects should play 
an important role and would avoid the breakdown of predictability, i.e., the disappearance 
of physical laws. The details of how this should actually happen is another mystery and 
probably different quantum theories of gravity would lead to different mechanisms for re- 
moving the singularities. 

If we accept that the idea of gravitation as a geometric phenomenon still persists at the 
quantum level ^ , with perhaps quantized areas and volumes in the fashion of loop quantum 
gravity [2| , it seems reasonable to expect that the quantum corrected gravitational dynamics 
could be described in terms of some effective action incorporating a number of regulating 
parameters (while keeping the matter sector untouched). In the absence of fully understood 



^ Note that string theory predicts the existence of fields with couplings that violate the equivalence principle 
[l| . For this reason, the idea of gravitation as a purely geometric phenomenon is explicitly broken in that 
context. 



quantum theories of gravity and of their corresponding effective actions, it would be desir- 
able to have a working model which, as an intermediate step between classical GR and the 
final quantum theory of gravity, could capture, at least qualitatively, some aspects of the 
sought non-singular theory of gravity. Stated differently, can we find a regulated gravita- 
tional theory (free from singularities) and as successful as GR at low energies? Obviously, 
such a theory would be very welcome from a phenomenological point of view and could 
provide new insights on fundamental properties of the geometry at very high energies. 

In this work we elaborate in this direction and propose a family of modified Lagrangians 
which departs from GR by quadratic curvature corrections 

f{R, R,,Rn = R + a^ + b^^ , (1.1) 

Kp rip 

where Rp ~ l~p^ is the Planck curvature, and show that for a wide range of parameters a and 
h they lead to non-singular cosmologies both in isotropic and anisotropic Bianchi I universes 
for all reasonable sources of matter and energy. In particular, we find that radiation domi- 
nated universes are always non-singular. The novelty of our approach, obviously, is not the 
particular Lagrangian considered, which is well known and naturally arises in perturbative 
approaches to quantum gravity. The new ingredient that makes our model so successful in 
removing cosmological singularities is the fact that we follow a first order (Palatini) formula- 
tion of the theory, in which metric and connection are assumed to be independent fields. In 
this approach, the metric satisfies second-order partial differential equations, like in GR, and 
the independent connection does not introduce any additional dynamical degrees of freedom 
(like in the Palatini version of GR [3i]). In fact, the connection can be expressed in terms of 
the metric, its first derivatives, and functions of the matter fields and their first derivatives. 
As a result, the theory is identical to GR in vacuum but exhibits different dynamics when 
matter and/or radiation are present. For the model 01. ip . this means that the dynamics is 
identical to that of GR at low curvatures but departures arise at high energies/curvatures. 
Since the equations are of second-order, there can not be more solutions in this theory than 
there are in GR. Therefore, the modified solutions that we find represent deformations (at 
the Planck scale) of the solutions corresponding to GR. Such deformations, as we will see, 
are able to avoid the big bang singularity by means of a bounce from an initially contracting 
phase to the current expanding universe. 

Our approach is motivated by previous studies on non-singular bouncing cosmologies 
initiated in [^ and continued in |5| and [6| . In j^ it was shown that the effective dynamics of 
loop quantum cosmology [7|, which describes an isotropic bouncing universe, can be exactly 
derived from an f{R) action with high curvature corrections in Palatini formalism. Different 
attempts to find effective actions for those equations followed that work but either failed [8| 
or are limited to the low-energy, perturbative regime |9|. The existence and characterization 
of bouncing cosmologies in the f{R) Palatini framework was studied in |5|, and the evolution 



of cosmological perturbations has been recently considered in [10|. The field equations of 
extended Palatini theories f{R,Q), in which the gravity Lagrangian is also a function of 
the squared Ricci tensor Q = R^^R^^^ , were investigated in p], where it was found that 
in models of the form f{R,Q) = f{R) + Q/Rp, the scalar Q is generically bounded from 
above irrespective of the symmetries of the theory. Unlike in the more conventional metric 
formalism. Palatini Lagrangians of the form f{R, Q) lead to second-order equations for the 
metric and, therefore, are free from ghosts and other instabilities for arbitrary values of the 
parameters a and b. 



The successful results of [^ and |5[ motivate and force us to explore scenarios with 
less symmetry to see if Palatini theories are generically an appropriate framework for the 
construction of non-singular theories. We will see that f{R) models which lead to bouncing 
cosmologies in the isotropic case also lead to anisotropic Bianchi I universes with expansion 
and energy density bounded from above. As we show here, however, such models have an 
unavoidable shear divergence, which occurs when the condition 9r/(-R) = is met. This 
important result implies that Palatini f{R) theories do not have the necessary ingredients 
to allow for a fully successful regulated theory in the sense defined above. Such limitation, 
however, is not present in f{R, Q) Palatini theories. We explicitly show that for the model 
(11. ip there exist bouncing solutions for which the expansion, energy density, and shear are 
all bounded. This model, therefore, avoids the well known problems of anisotropic universes 
in GR, where anisotropics grow faster than the energy density during the contraction phase 
leading to a singularity, which can only be avoided by means of matter sources with equation 
of state w = P/p> 1 [ll|. 

The content of the paper is organized as follows. In section |Tl] we summarize the field 
equations of Palatini f{R, Q) theories with a perfect fluid, which where first derived and 
discussed in |6|] . In section IIIII we obtain expressions for the expansion and shear in both 
f{R) and f{R,Q) theories. Section |IV] is devoted to the analysis of f{R) theories in 
isotropic and anisotropic scenarios, paying special attention to the possible existence of 
isotropic bouncing solutions which are not of the type dnf = 0. In section |V] we study the 
model (11. ip and characterize the different bouncing solutions according to the values of the 
Lagrangian parameters a and b, and the equation of state w. We end with a brief discussion 
and conclusions. 



II. FIELD EQUATIONS 

The field equations corresponding to the Lagrangian (11.11) can be derived from the action 

S = ^J d'xy^fiR, Q) + SU9a^, ^) (2.1) 

where R ^ g'^^R,,, Q ^ R,,R^^, R,, ^ -d,Tl + d,T'^^ + T'^Xpx ' ^tp^U' ^U is the inde- 
pendent connection, and \& represents generically the matter fields, which are not coupled 
to the independent connection. Variation of the action with respect to the metric leads to 

f 
faRpu — ■^fi'/ii^ + '^fqRfj.aR^u = /t Tf^u , (2.2) 

where fji = d^f and /g = dqf. Variation with respect to the independent connection gives 

V/. [V^ ifng^'^ + 2fQR^n] = (2.3) 



For details on how to obtain these equations see [6|. The connection equation (12. 3p can be 
solved in general assuming the existence of an auxiliary metric hai3 such that (12. 3 p takes the 
form V/3 [\/—hh^'^~\ =0. If a solution to this equation exists, then r(^^ can be written as 
the Levi-Civita connection of the metric /i^jy. When the matter sources are represented by 



a perfect fluid, T^i, = (p + P)u^Ui, + Pg^vi one can show that h^^ and its inverse h^" are 
given by 



h 
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In terms of /i^i, and the above deflnitions, the metric field equations (12. 2p take the following 

form 
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RAh) = ^ 
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In this expression, the functions /, Ai, and A2 are functions of the density p and pressure 
P. In particular, for our quadratic model one finds that R = H?{p — 3P) and Q = Q{p, P) 
is given by 



hQ 
2Rp 



f , R 



P f2 



k'p + ^ + -^a + 

2 8b ■'^' 
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<f^/«)-\/(f^/«)-^^^^ 



(2.12) 

where f = R + aR^/Rp, and the minus sign in front of the square root has been chosen to 
recover the correct limit at low curvatures. 



In what follows, we will use (12. lip to find equations governing the evolution of physical 
magnitudes such as the expansion, shear, matter /energy density, and so on. Note that (12. lip 
is written in terms of the auxiliary metric /i^i,, not in terms of the physical metric g^y. In 
terms of g^i,, eq. (12.111) would be much less transparent and more difficult to handle. As 
we will see in the next section, working with (12. lip will simplify many manipulations and 
will allow us to obtain a considerable number of analytical expressions for all the physical 
magnitudes of interest. 



III. EXPANSION AND SHEAR 



In this section we derive the equations for the evolution of the expansion and shear for 
an arbitrary Palatini f{R,Q) theory. We also particularize our results to the case of f{R) 



theories, i.e., no dependence on Q. We consider a Bianchi I spacetime with physical hne 
element of the form 

3 



ds' 



g^ydx^dx" = -dt^ + ^a^(t)((ix* 



(3.1) 



4 = 1 



In terms of this line element, the non-zero components of the auxiliary metric h^y are the 
following 
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The relevant Christoffel symbols associated with h^u are the following 
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The non-zero components of the corresponding Ricci tensor are 






(3.7) 
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where if^ = d^/O'k- These expressions define the Ricci tensor R^y{h) on the left hand side of 
eq. (12. lip . For completeness, we give an expression for the corresponding scalar curvature 



R{h) = ^ 



2Y,H, + Y,Hl+[Y,Hu 



Vt 



n 



's n 



Y.H. 



Q 3QS 
'q'2Q^ 



(3.9) 

From the above formulas, one can readily find the corresponding ones in the isotropic, 
fiat configuration by just replacing Hi — )• H. For the spatially non-flat case, the Rtt{h) 
component is the same as in the flat case. The Rij{h) component, however, picks up a new 
piece, 2K'jij, where jij represents the non-flat spatial metric oi gij = af'jij- The Ricci scalar 
then becomes R{h) -^ R^=%h) + g. 



A. Shear 



From the previous formulas and the field equation fl2.1ip . we find that i?/ — Rj^ = (no 
summation over indices) leads to 



Ri' — Rj^ = — 



Hij + Hij <^ 2_, ^fc + 7 



_ , fi S 

k 



, (3.10) 



where we have defined Hij = Hi — Hj. Using the matter conservation equation for a fiuid 
with constant equation of state P = wp, 



P 



[l + w)pJ2Hk, (3.11) 



the above equation can be readily integrated (for this reason we consider constant equations 
of state throughout the rest of the paper). This leads to 

1 1 1 

5'2p(i+iu) pTTT^ 

Hij = Cij —3 = Cij- — (3.12) 

iU Ai — A2 

where the constants Cij = —Cji satisfy the relation C12 + C23 + C31 = 0. It is worth noting 
that writing explicitly the three equations (I3.12p and combining them in pairs, one can write 
the individual Hubble rates as follows 

1 

6 p (! + «') 

Hi = - + (C12 — C31) 



3 ' '^ -^^^Ai-A: 



2 



1 



9 n(i+™) 

H2 = - + {C23 - C12) f^^j^^ (3.13) 

1 
9 p(i+"') 

H3 = - ~^ (C*31 " C23) 



3 ' -" ^' Ai - A 



2 



where 9 is the expansion of a congruence of co moving observers and is defined as ^ = ^^ H^. 
Using these relations, the shear a"^ = J2i {Hi — |) of the congruence takes the form 



,_ p^ ic!2 + ch + cL 



^2 r \^i2 ' —-^d ' ^6iJ In 1 /i\ 

- = (aTTa^i 3 • P-"' 



where we have used the relation (C12 + C23 + C31) = 0. 



B. Expansion 

We now derive an equation for the evolution of the expansion with time and a relation 
between expansion and shear. From previous results, one finds that 

k \ k / k \ / 



In terms of the expansion and shear, this equation becomes 



Cti^-y + y (^1+2^1) , (3.16) 

where we have defined 

Ai = -(1 + ^)P^ (3.17) 

The right hand side of fl3.16p is given by G^j,^ = r^^, — \h^yh°^I^Taji, being t^^ the right hand 
side of (12.111) . A bit of algebra leads to the following relation between the expansion, shear, 
and the matter 



Note that once a particular Lagrangian is specified, an equation of state P = wp is given, 
and the anisotropy constants Cij are chosen, the right hand side of Eqs. fl3.14p and fl3.18p 
can be parametrized in terms of p. This, in turn, allows us to parametrize the Hi functions 
of fl3.13p in terms of p as well. This will be very useful later for our discussion of particular 
models. 

In the isotropic case (a^ = 0,6* = 3a/ a = 371) with non-zero spatial curvature, fl3.18p 
takes the following form 

1 \f + K\p + 3P)-^^] , , 

n^ = —- — —^ \ "' ^ 3.19 

6(Ai-A2) [l + fAi]' ^ ' 

The evolution equation for the expansion can be obtained by noting that the Rij equa- 
tions, which are of the form Rij = {Q/2S)gij [...] = (//2 + H?P)gij/ Ki, can be summed up 
to give 

Using the relations tl = —{l + w)pQp6, S = —{l + w)pSp6, and Q = {l + w)'^p6'^[Qp + pQpp] — 
(1 + w)pQp9, the above expression turns into 



[2+3Ai]^+ 



2 + (2 - 3w)A, + 3A2 - (1 + ^)p ( 1 + ^Ai j (^ - y 



e 



2 3 [/ + 2k^P] 



A1-A2 ' 
(3.21) 

where we have used the definition (13 . 1 7p and have defined the quantity 

A2^(1+«;)V^ (3.22) 

Note that the function 6 can also be plotted as a function of p. In the isotropic, non-fiat 
case the evolution equation for the expansion {6 = 3H) can be obtained from (I3.2ip by just 
replacing the term [/ + 2h?P] on the right hand side by [/ + 2n^P — AKKi/a?]. 



C. Limit to f{R) 



We now consider the limit /g — t- 0, namely, the case in which the Lagrangian only depends 
on the Ricci scalar R. Doing this we will obtain the corresponding equations for shear and 
expansion in the f{R) case without the need of extra work. This limit can be obtained from 
eqs.f l2.6p to (12. 9p by taking /q — )■ in those definitions. One then finds that 



Ai ^ /i? , A2 -> 

S ^ n^fR 



Equation (I3.12p turns into 



which leads to 



1 






,(l+m) 



Ir 



(3.23) 
(3.24) 



(3.25) 
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(3.26) 
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The shear in thus given by 
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2_p^{ci, + ci, + a 



3lJ 



pR 3 

where C12 + C23 + C31 = 0. The relation between expansion and shear now becomes 



(3.27) 



— I + -A1 
3 V 2 ' 



f + K'^p + W) a^ 
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(3.2^ 



where Ai is given by (I3.17P but with Vt replaced by Jr. In the isotropic case with non-zero 
K we find 



W 



6/. 



R 



1 + |A, 



The evolution equation for the expansion is now given by 

_ 3 [/ + 2«:2p - 1^] 



[2 + 3Ai]e + 2 + (2 - 3w)Ai + 3A2 



/ 



R 



(3.30) 



where A2 is defined as in (13.220 but with Vt replaced by Jr. 
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IV. ISOTROPIC AND ANISOTROPIC BOUNCING f{R) COSMOLOGIES 

An isotropic and homogeneous cosmological model experiences a bounce when the 
Hubble function "H^ vanishes (see FigH]), thus defining a minimum of the expansion factor 
(see Figj2]). According to the formulas derived in previous sections, isotropic bouncing 
f{R) cosmologies occur either when the denominator of (I3.29P blows up to infinity or 
when the numerator [/ + K^(p + 3P) — ^^7^] vanishes. The divergence of the denominator 
only depends on the form of the Lagrangian, whereas the vanishing of the numerator also 
depends on the value of the spatial curvature K. In order to characterize the anisotropic 
bouncing models, it is convenient to study first the isotropic case. For this reason, we will 
focus first on the existence of divergences in the denominator and will postpone until the 
end the other case. 



A. Divergences of Ai and importance of anisotropies. 

As can be easily verified from the definition of Ai in (I3.17p . the existence of divergences in 
the denominator of ( 13.290 can only be due to the vanishing of the combination /ij(-R/R_R— /h): 

^^ ^ il + w)il-3w)K'pfnR ^^^^^ 

fR{RfRR — Ir) 

The Lagrangian that reproduces the dynamics of loop quantum cosmology with a massless 



R 
12Rc 



scalar, which is well approximated by the function f{R) = — J dRtanh j -^ In 

satisfies the condition //j = at i? = 12Rc, where Re is a scale related with the Planck cur- 
vature Rp, thus leading to a divergence of Ai at that point. One can construct other models 
with simple functions such as f{R) = R + aR^/Rp or f{R) =R + R^/Rp{a + h\n[R^/Rl]) 
which also have bounces when fp = 0. The fp = bouncing condition seems to be 
quite generic and arises even when one tries to find models which satisfy the condition 
R/rr — /r = at some point. An illustrative example is the model f{R) = Rp{e^^^^ — 1), 
which leads to fp = e^/^^ and Rfpp-fp = e^/^p{R-Rp)/Rp, which vanishes at R = Rp. 
In this model one either finds a divergent "H^, due to the vanishing of the denominator 
of 03.29P for w < 1/3, or a bounce when the density approaches the limiting value 
k^Pb = 2Rp/{3w — 1) for w > 1/3. This bounce occurs as R/ Rp ~ ln[l — p/ pp] — ^ —00, 
which corresponds to /^j — )■ and, therefore, lies in the standard class of bouncing models. 

The importance of finding f{R) models for which the bounce occurs when fp 7^ becomes 
apparent when one studies anisotropic (homogeneous) scenarios. In these cases, the shear 
diverges as ~ 1//^, as is evident from (13.270 . This shows that any isotropic bouncing 
cosmology of the /r = type will develop divergences when anisotropies are present. And 
this is so regardless of how small the anisotropies are initially. It is worth noting that 
eventhough cr^ diverges at /i^ = 0, the expansion and its time derivative are smooth and 
finite functions at that point if the density and curvature are finite. In fact, from (I3.28P and 
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Figure 1: Representation of the Hubble function in terms of p for the model /(i?) = R — B? /2Rp 
and u; = for K < 0, K = 0, and K > Q. 
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Figure 2: Time evolution of the expansion factor for the model f{R) = R — R? /2Rp and w = Q 
for X > 0, -fC = 0, and K <{) (solid curves from left to right). From left to right, we see that the 
universe is initially contracting, reaches a minimum, and then bounces into an expanding phase. 
The dashed lines visible near the bounces represent the expanding solutions of GR, which begin 
with a big bang singularity and quickly tend to the non-singular solutions. 



fl3:30|l we find that^ 



01 



2(^*12 + ^"23 + C'si; 



-Rq 



(2-3^) 2 

° ' 2(1 -3w) 



-el + 



(1 + w)(l — 3w)k'^Po 



Po 



2 

l-l-tu 



(4.2) 
(4.3) 



^ Note that the case w — 1/3 must be excluded from the analysis because in that case the theory behaves 
like GR with an effective cosmological constant and the manipulations that lead to (|4.2p and (|4.3p are not 
valid. 
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0.10 
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Figure 3: Evolution of the expansion 9^ with the amount of anisotropy. The soHd (blue) line 
represents the isotropic case. As the anisotropy grows, the first local maximum of 9^ disappears 
(dashed curves) . This indicates that anisotropies can overcome the repulsive forces of the modified 
f{R) Lagrangian. 



where the subindex denotes the point at which Jr = (where the shear diverges). It is 
worth noting that in GR ^ < always, whereas in f{R) the point /r = is characterized by 
(14. 3p . which may be positive, negative, or zero. If the anisotropy is sufficiently small, which 
is measured by the constant (C'^2 + C'23 + C'|i) in (14. 2p . then 6 may be positive. This indicates 
that some repulsive force is trying to halt the contraction. However, if the anisotropy is too 
large, then it can dominate the expansion and keep ^ < at all times (see Figj3] and note 
how the first local maximum tends to disappear in the upper curves as the anisotropy grows). 

In Figj3] we find that there exist anisotropic solutions for which ^ = at densities 
beyond the point /^j = 0, which sets the bounce of the isotropic case. One could 
thus be tempted to claim that for universes with low degree of anisotropy bouncing 
solutions really exist if we allow for slightly negative values of fn, which lead to ^ = 0. 
However, the shear divergences of these anisotropic models at //? = are physically 
unacceptable because any detector crossing the singularity would be ripped apart by 
the infinite tidal forces (see |l2| for a nice discussion on divergences and singularities in 
cosmology). Moreover, from Eqs. (I3.26P it is easy to see that the Kretschman scalar 
R^.uapR^"'"' = ^J2iiHi + HfY + HfHl + HfHl + HlHl) diverges at least as ~ 1//^, which 
is a clear geometrical pathology. Additionally, the vanishing of f^ suggests that the field 
equations may not be valid for negative values of fa because the conformal transformation 
needed to solve for the connection becomes ill-defined at /r = 0, which seems to be a 
generic problem of anisotropic models in modified theories of gravity [l3|. Note also that 
the evolution of inhomogeneous perturbations in isotropic models develops divergences 
when fpt vanishes 10 . 



When the bounce is due to the vanishing of RfuR — j'r aX R = Rb (with j'r^ Q aX that 
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point), then the shear is finite and the expansion is given by 

Ol ~ (RfRR - IrT ^ (4.4) 

7 + n\l + 'iw)pB (J^ 



k _ ^RsfRRR 

^ ~ 3{1 + 3){1-3w)k'^Pb 



Vr ^ 2 



(4.5) 

Rb 



The fact that ^^ = at Rb imphes that the density reaches a maximum at that point 
(recall the conservation equation: p = —6{p + P)). Also, since in this case the shear is 
finite, this family of bouncing f{R) models seems to be the right family of Lagrangians to 
construct non-singular models. However, as we show next, there are no Lagrangians of this 
type able to recover GR at low curvatures. 



1. Non-existence of RfRR — fR = models. 

The existence of bounces in the isotropic case is due to the unbounded growth of 
i 1 + I Ai j . One may try to build bouncing models by defining an always positive func- 
tion g{R) which has a divergence at i? = Rp such that 

f ^x\ {fRR\Q{^ + w)f-(l + 3w)RfR]-fl) ,^ , 

V 2 y jR[RjRR - JR) 

Given the function g{R), one can find the Lagrangian f{R) that generates the corresponding 
bouncing Universe by just solving a second order differential equation. The function g{R) 
also needs to satisfy the condition g{R) ~ 1 as -R — ?■ to force f{R) ~ -R in that limit. 
Simple manipulations of (14. 6 p lead to 

Irr ^ [2 - gjR)] .4 . 

/2 Q{i + w)f-[l + 3w + g{R)]RfR- ^-' 

Since in GR R > ii w < 1/3 and i? < if w > 1/3, we may perform the change of variable 
f{R) = ±Roe^^^\ which leads to /r = Xr/, Jrr = (Xrr + A|j)/, and allows us to rewrite 
(gTj) as follows 

^RR + A| ^ [2 - g{R)] .4 g. 

A| Q{1 + w) ~ [1 + 3w + g{R)]R\R ' ^'^ 

By construction, the function g{R) goes like g{R) ~ 1 at low curvatures, then may change in 
an unspecified way though always being positive at intermediate curvatures, and finally blows 
up to infinity at i? = Rp, which sets the high- curvature scale. Since g{R) grows unboundedly 
near Rp, we see that the denominator of (14.81) could vanish at some point. This is in fact what 
one finds systematically when using a numerical trial and error scheme to find f{R) bouncing 
models. We now show that this always occurs for any function g{R) satisfying the conditions 
required above. Since at low curvatures we demand f{R) ~ R, which implies Xr ~ 1/-R, it 
follows that the denominator of (14. 8 p is Den ~ 6(1 +w) — {1 + 3w + 1)1 = (4 -|- 3w), which 
is positive for all reasonable matter sources {w > —4/3). After this initial positive value, 
since g{R) > will grow as RXr remains positive^, unavoidably we will have Den = at 
some later point. Then: 



^ Note that the product RXr is initially positive and can only change sign if A^ vanishes at some point, 
which would force fa — at that point. 
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• If Ar 7^ when Den = 0, then g{R) and Xr are finite whereas A/j^ — ;• oo at that point. 
However, since g{R) is finite, the divergence of Xrh cannot imply a cosmic bounce, 
since by construction that only happens when g{R) diverges. Therefore, this case does 
not correspond to a bounce. 

• If we admit that g{R) can indeed go to infinity, it follows that that must be the only 
point at which Den = 0. This requires that the product g{R)XR be finite at Rp, which 
implies that A/j — )■ at Rp to exactly compensate the divergence of g{R) and give a 
final result which exactly cancels with the 6(1 + w) of the denominator of (14. 8p . Note 
that in this case the left hand side of (14. 8 p diverges as l/X'p and the right hand side 
goes like —g{R)/ zero. 

This shows that the bouncing condition g{R) — )■ oo at Rp can only be satisfied if Xr vanishes 
at that point, which implies that fp = and excludes the possibility of having RfpR — fR = 
as the condition for the bounce. 



B. Vanishing of the numerator of Ti^ 

In the previous subsection, we concluded that the denominator of "H^ can only diverge if 
fp = 0. We now investigate if there exist some other mechanism able to generate isotropic 
bouncing models. We begin by noting that the bounce should occur when f^+{l+3w)K^pB — 
QKf^/a\ = 0. Using the well-known relation 

Rfp - 2/ = k'T , (4.9) 



which follows from the trace of the field equations of Palatini f{R) theories, we find that 

^B - ^ j /I = -3(1 + w)k^pb (4.10) 



a. 



Now, since at low curvatures //? ~ 1 > and must remain positive always (to avoid a bounce 
of the type fp = 0), at the bounce we have {R^ — 12K/a^) < for all w > —1. Since at 
low densities i? ~ (1 — 3w)k'^p is positive for w < 1/3, the negative sign of {R^ — 12K/a%) 
implies that for K < the rate of growth of R with p must vanish and change sign at some 
point before the bounce. Using eq. (l4.9p . we find that 

JR — J^JRR 

A change of sign in dpR implies a divergence in the denominator of this last equation, which 
means that fp^oo and/or fpp — )■ — oo. In none of those cases the theory is well defined 
beyond the divergence, which implies that R is monotonic with p. Therefore, for i^ < 
the only hope is a bouncing model with w > 1/3 because for such equations of state R < 
always. For i^ > this constraint can, in principle, be avoided. 

Let us now focus on the case K = 0. We can parallel the strategy followed in the 
previous section and build f{R) models starting with a function g{R) which goes like R at 
low curvatures and has a zero at i? = -Rp such that 

^2 ^ 9{R) ^^ _ (4_^2) 



3(1-3^)/r(i + ^) 
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The function g{R) determines a first order differential equation, 2g{R) = (1 + 3w)Rfji — 
3(1 + w)f, from which f{R) can be easily obtained as 

/(i?) = --^^ r dx^ (4.13) 

^^ ' {l + 3w) J xi+^ ^ ' 

where 7 = | g^ . Though this is a convenient method for model building, a trial an error 
analysis does not lead to any successful model^. Numerically, we find that either an //? = 
bounce occurs or that the denominator of "H^ vanishes before the zeros of g{R) can be 
reached, which leads to a singularity. 

When K j^ 0, the above method can also be applied, though the resulting differential 
equation becomes highly non-linear and the solutions can only be found numerically. The 
results are similar to the case K = 0. We systematically find that the models with a hope 
to lead to a bounce are those for which /« — )■ at some point. As a result, the spatial 
curvature term —QKfji/a? is suppressed in that region and becomes negligible, giving rise 
to a bounce of the type f^ = 0. Though a rigorous proof similar to that given in the case 
of RJrr — Jr = models is not yet available, we believe that no models of this type which 
recover GR at low curvatures exist. 



C. Conclusions for f{R) models 
Using Eqs. 03.28P and 03.291) . the expansion can be written as follows 

02 = gH-' + ^ i^, (4.14) 

where H represents the Hubble function in the K = isotropic case, and cr^ is defined in 
(I3.27p . From this representation of the expansion, it is clear that the only way to get a 
true bouncing model without singularities is by satisfying the condition R/rr — /r = 0, 
which would generate a finite shear, a divergent denominator in the second term of (14.141) . 
and hence a vanishing expansion. However, we have explicitly shown that such condition 
can never be satisfied. Moreover, even if the numerator of H^ could vanish and produce a 
different kind of isotropic bouncing models, in the anisotropic case the expansion would not 
vanish and, therefore, that could not be regarded as an anisotropic bounce. For all these 
reasons, it follows that Palatini f{R) models do not have the necessary ingredients to build 
a complete alternative to GR free from cosmic singularities. 

V. NONSINGULAR UNIVERSES IN f{R,Q) 

The previous section represents a no go theorem^ for the existence of non-singular Palatini 
f{R) models able to produce a complete alternative to GR in scenarios with singularities. 



^ Among many others, we considered families of models characterized by functions such as g{R) = i?(l 
R^/R'pY' and g{R) = R{1 - {R/Rpf lni?V^p)- 
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This is so at least for universes filled with a single perfect fluid with constant equation of state. The 



consideration of fluids with varying equation of state lj| or with anisotropic stresses, see for instance 
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Though the isotropic case greatly improves the situation with respect to GR, the anisotropic 
shear divergences kill any hopes deposited on this kind of Lagrangians. The most natural 
next step is to study the behavior in anisotropic scenarios of some simple generalization of 
the f{R) family to see if the situation improves. Using the Lagrangian fll.ip . we will show 
next that completely regular bouncing solutions exist for both isotropic and anisotropic 
homogeneous cosmologies. 



A. Isotropic Universe 

Consider Eq. fl3.19l) together with the definitions fl2.6p - fl2.9l) particularized to the f{R, Q) 
Lagrangian (11.11) . In this theory, we found that R = K^{p — SP) and Q = Q{p, P) is given 
by (I2.12p . From now on we assume that the parameter h of the Lagrangian is positive and 
has been absorbed into a redefinition of Rp, which is assumed positive. This restriction is 
necessary (though not sufficient) if one wants the scalar Q to be bounded for w > —1. The 
Lagrangian then becomes /(-R, Q) = R + aR^/Rp + Q/Rp. When b/Rp > 0, positivity of 
the square root of eq. f l2TT2|) establishes that there may exist a maximum for the combination 
p + P. 

The first difficulty that we find is the choice of sign in front of the square root of Eq. (l2.8p . 
In order to recover the f{R) limit and GR at low curvatures, we must take the minus sign. 
However, when considering particular models, which are characterized by the constant a and 
an equation of state w, one realizes that the positive sign and the negative sign expressions 
for A2 may coincide at some high curvature scale, when the argument of the square root 
^yX^ — k'^{p + P) vanishes. When this happens, one must make sure that the function A2 
at higher energies is continuous and differentiable. These two conditions force us to switch 
at that point from the negative to the positive sign expression (see FigH] for an illustration 
of this problem), which then defines a continuous and differentiable function on the physical 
domain. Bearing in mind this subtlety, one can then proceed to represent the Hubble func- 
tion for different choices of parameters to determine whether bouncing solutions exist or not. 

We observe that for every value of the parameter a there exist an infinite number of 
bouncing solutions, which depend on the particular equation of state w. The bouncing 
solutions can be divided in two large classes: 

• Class I: a > 0. 

In this case, the bounce occurs at the maximum value reachable by the scalar Q. This 
happens when the argument of the square root of (I2.12p vanishes. In general, the 
density at that point satisfies 

^Vq^.. _ 1 + 5w - 2a(l - 3w) - V8(l + w){2w - a(l - 3w)) 

Rp (l + 2a)2(l-3u;)2 ' ^ ' ' 

The bounce occurs at that density for all equations of state satisfying the condition 

w > ^_ = ^ (5.2) 



15|, could affect the dynamics adding new bouncing mechanisms, and potentiahy restrict the range of 



applicability of this conclusion. 
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w=-l/3 
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Figure 4: Illustration of the need to combine the two branches of A2 to obtain a continuous and 
differentiable curve. The branch that starts at the origin has the minus sign in front of the square 
root (solid line). When the square root vanishes, the function must be continued through the 
growing dashed branch, which corresponds to the positive sign in front of the square root. The 
matching point of the radiation universe {w = 1/3) occurs at K?p = Rp/6 and appears highlighted 
in green (online only). 

which follows from the requirement of positivity of the argument of the square root 
of (15.11) . Note that, except for a = 0, the case w = Wmin is not contained in the 
set of bouncing solutions. From (I5.2p it follows that a radiation dominated universe, 
w = 1/3, always bounces for any a > 0. In fact, when w = 1/3, we find that (15. ip 
must be replaced by 

-,m=i/3 = ^ ^ ■ (5-3) 



K^p, 
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Note that this last expression is independent of the value of a and, thus, holds also 
for the case a < 0. This was to be expected since the coefficient a multiplies the 
quadratic term R^/Rp, which is zero in a radiation dominated universe. Remarkably, 
this implies that all radiation dominated universes in the family of Lagrangians 
considered here always lead to a big bounce. This clearly demonstrates that f{R,Q) 
theories posses interesting dynamical properties that cannot be reproduced by any 
f{R) Palatini Lagrangian [l6|. The modified dynamics in the f{R) case is generated 
by new terms that depend on the trace T = — (1 — 3w)p, which do not produce any 
effect in a radiation scenario. 



Class II: a < 0. 

This case is more involved and must be divided in several intervals. In general, the 

bounce occurs at a density given by the following expression 



k^Pb 
Rp 


( l+ew-2a(l-3w)-3-s/w(2+3w)-a{l+w){l-3w) 


) (l+a)(l+4a)(l-3«>)2 

1 «^PQma^ 



Rp 



ii w < Wo 
if w > Wq 



(5.4) 
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where Wo represents the value of w at which the two branches of k'^pb/Rp coincide. 
The generic expression for wq as a function of a is very comphcated, though its 
computation for a given a is straightforward. Note that the curve defined by Eq.f l5.4p 
is smooth and differentiable with respect to w even at Wq. It is important to note that 
Wq is always negative. This means that the bouncing solutions that occur at pq^^^ 
can be extended to negative values of w until the value wq. As of that point, the range 
of bouncing solutions is extended to even more negative values of w through the new 
branch w < wq of Eq.f l5.4p . What happens before and after wq to make that particular 
equation of state so relevant? The answer is as follows. For w > Wq, the bounce 
occurs at a density for which Q is maximum (when the square root of fl2.12p vanishes). 
For w < Wo, the bounce occurs at a density for which the function Ai — A2 vanishes. 
At w = Wq we find that Q reaches its maximum at the same density as Ai— A2 vanishes. 

How far into the negative axis can w be extended beyond the matching point wq? The 
answer depends on the value of a. We split the a < axis in five elements: 

- Case Ila: -1/4 < a < 

The values of w in this interval are restricted by the argument of the square root 
of f l5.4p for w < Wq. We thus find that 



We see that when a = we find agreement with the discussion of Case 1. As a 
approaches the limiting value —1/4, the bouncing solutions extend up to tt; — )• 
— 1/3. However, since the branch w < Wq oi (15.41) is singular at a = —1/4, that 
particular model must be studied separately. 

Case lib: a = -1/4 

In this case, the density at the bounce is given by the following expression 

Rp I ^tm^ if «; > -| ^ ' 

which is always finite except for the limiting value w = —1/3. Thus, bouncing 
solutions exist for any w within the interval —1/3 < w < 00. 

Case lie: -1/3 < a < -1/4 

Though in this interval the argument of the square root in (15. 4p is always positive, 

we observe numerically that the bouncing solutions cannot be extended beyond 

the value w < —1, where pb reaches a maximum. Therefore, in this interval 

we find that the bouncing solutions occur if —1 < w < 00, where w = — 1 is 

excluded. 

Case lid: -1 < a < -1/3 

Here we also find that the negative values of w cannot be extended beyond 
w < —1. Surprisingly, we also find restrictions for w > 1 which are due to the 
existence of zeros in the denominator of H^. Due to the algebraic complexity of 
the functions involved, it is not straightforward to find a clean way to characterize 
the origin of those zeros. However, numerically we find that they arise when 
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w ~ (a + /3a)/(l + 3a)^, where a = 1.1335 and (3 = —3.3608 (this fit is very good 
near a ^ —1/3 and shghtly worsens as we approach a = —1). Summarizing, the 
bouncing solutions are restricted to the interval —l<w< (a; + /3a)/(l + 3a)^ > 1. 
This expression agrees in the limiting value a = —1/3 with the expected values 

— 1 < w < oo of the case lie. The case a = —1 must be treated separately, 
though it does not present any undesired feature. 

Note that in this interval one finds the case a = —1/2, which is singular according 
to (15. ip and must be treated separately. We find that Eq. (15.11) must be replaced 
by k'^Pq^^^ = 1/(4 + 4:w). Other than that, this case satisfies the same rules as 
the other models in this interval. 

— Case lie: a < — 1 

Similarly as the family a > 0, this set of models also allows for a simple 
characterization of the bouncing solutions, which correspond to the interval 

— 1 < w < a/{2 + 3a). In the limiting case a = —1 we obtain the con- 
dition —l<w<l (compare this with the numerical fit above, which gives 

— l<u'<1.12). In that case, the density at the bounce is given by 

Rp - |!lig^ if w;>-| ^^■'> 

For a < — 1, the equations of state that generate bouncing solutions get reduced 
from the right and approach —l<w< 1/3 as a — )■ — oo, with the case w = 1/3 
always included. 

B. Anisotropic Universe 

Using Eqs. (I3.18P and (I3.19p . the expansion can be written as follows 

where H represents the Hubble function in the K = {) isotropic case. To better understand 
the behavior of 6'^, let us consider when and why H^ vanishes. Using the results of the 
previous subsection, we know that H^ vanishes either when the density reaches the value 
pQmax or when the function Ai — A2 vanishes. These two conditions imply a divergence in 
the quantity (1 + |Ai)^, which appears in the denominator of H^ and, therefore, force the 
vanishing of H^ (isotropic bounce). Technically, these two types of divergences can be easily 
characterized. From the definition of Ai in (I3.17p . one can see that Ai ~ dpQ/Q. Since 
n = a/Ai(Ai — A2), it is clear that Ai diverges when Ai — A2 = 0. The divergence due to 
reaching Pq„„^ is a bit more elaborate. One must note that dpQ ~ dpAi ~ dpA2 and that 
dpAi^2 contain terms that are finite plus a term of the form dpX, with A given by (12. 9p . In 
this A there is a Q— term hidden in the function f{R, Q), which implies that dpX ~ dpQ/Rp 
plus other finite terms. From the definition of Q it follows that dpQ has finite contributions 
plus the term dp^/y/^, where $ = (1 + (1 + 2a)R/Rp)'^ - Ak^{p + P)/Rp, which diverges 
when $ vanishes. This divergence of dpQ indicates that Q cannot be extended beyond the 
maximum value Qmax- 

Now, since the shear goes like a^ ~ l/(Ai — A2)^ [see Eq. (l3.14p ]. we see that the condition 
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Ai — A2 = implies a divergence on o"^ (though 9^ remains finite). This is exactly the same 
type of divergence that we already found in the f{R) models. In fact, the decomposition 
(15. 8 p is also valid in the /(-R) case, where A2 — )■ and Ai — )■ fji (see Eq. (l4.14p ). Since in 
those models the bounce can only occur when /^ = 0, which is equivalent to the condition 
Ai— A2 = 0, there is no way to achieve a completely regular bounce using an f{R) theory. On 
the contrary, since the quadratic /(-R, Q) model (II. ip allows for a second mechanism for the 
bounce, which takes place at Pq^^^ , there is a natural way out of the problem with the shear. 

When the density reaches the value Pq„^^, we found in the previous subsection that the 
combination Ai — A2 is always greater than zero except for the particular equation of state 
w = wq (recall that wo was defined as the matching condition in eq. (l5.4p . and represents the 
case in which pg^^^ is reached at the same time as Ai — A2 = 0). Therefore, for any w > Wq 
the shear will always be finite at pq^^^- Moreover, since at that point the denominator 
(1 + |Ai)^ blows up to infinity, it follows that the expansion vanishes there, which sets a 
true maximum for p like in the isotropic case. 

At this point one may wonder about the consequences of the divergence of dpQ at pq^^^ 
for the consistency of the theory. This question is pertinent because the connection that 
defines the Riemann tensor involves derivatives of Vt and hence of Q. In this sense, it 
should be noted that because of the spatial homogeneity only time derivatives of such 
quantities need be considered. We are thus interested in objects such as dtQ and higher 
time derivatives. One can check by direct computation that dtQ = {dpQ)p yields a finite 
result because the divergence of dpQ is exactly compensated by the vanishing of p, which 
is due to the vanishing of the expansion at the bounce. Explicit computation of higher 
derivatives of Q and other relevant objects (such as f2, S, ... needed to compute the 
components of the Ricci tensor) shows that all them are well behaved at the point of the 
bounce^. This guarantees that the bounce is a completely regular point that does not spoil 
the well-posedness of the time evolution nor the disformal transformation needed to relate 
the physical and the auxiliary metrics g^i, and hp^,, respectively. 

Summarizing, we conclude that the Lagrangian (ll.ip leads to completely regular bouncing 
solutions in the anisotropic case for w > ^^ ii a > 0, ior wq < w < 00 if —1/3 < a < 0, 
for Wq < w < {a + /3a)/(l + 3af if -1 < a < -1/3, and for -1/3 < w < a/{2 + 3a) if 
a < —1, where wq < is defined using (15. 4p and its corresponding subcases. These results 
imply that for a < the interval < w < 1/3 is always included in the family of bouncing 
solutions, which contain the dust and radiation cases. For a > 0, the radiation case is always 
non-singular too. 



C. An example: radiation universe. 

As an illustrative example, we consider here the particular case of a universe filled with 
radiation. Besides its obvious physical interest, this case leads to a number of algebraic 



^ This same reasoning can be used to confirm the pathological character of the other type of bounce, the 
one characterized by the condition Ai — A2 = 0. 
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simplifications that make more transparent the form of some basic definitions 



Q 



A 



Ai 



3Rp 




(5.9) 
(5.10) 
(5.11) 



It is easy to see that the coincidence of the two branches of A2 occurs at K^p = Rp/6. 
Therefore, the physical A2 must be defined as follows (see FigJl]) 
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This definition by parts unavoidably obscures the representation of other derived quanti- 
ties. Nonetheless, it is necessary to obtain continuous and differentiable expressions for the 
physical magnitudes of interest such as the expansion and shear (plotted in Figs. [5l [6], and 
CD. It is easy to see that at low densities (ES]) leads to Q k, 4(kV)V3 + ?>2(K^pf/9Rp + 
320{k'^ p)'^/ 27 Rp + . . ., which recovers the expected result for GR, namely, Q = 3P^ + p^. 
From this formula we also see that the maximum value of Q occurs at n'^Pmax = 3i?p/16 
and leads to Qmax = 3-Rp/16. At this point the shear also takes its maximum allowed value, 
namely, ct"^^^ = ^/3JTQRp (Cfg + C2^ + C|i), which is always finite. At p^ax the expansion 
vanishes producing a cosmic bounce regardless of the amount of anisotropy. 



VI. DISCUSSION AND CONCLUSIONS 

In this work we have shown that simple modifications of GR with high curvature 
corrections in Palatini formalism successfully avoid the big bang singularity in isotropic and 
anisotropic (Bianchi-I) homogeneous cosmologies giving rise to bouncing solutions. And this 
type of solutions seem to be the rule rather than the exception. The f{R, Q) model (II. ip in 
Palatini formalism is just an example. This type of models is motivated by the fact that the 
effective dynamics of loop quantum cosmology [TJ is described by second-order equations and 
by the need to go beyond the dynamics of Palatini f{R) theories, which cannot avoid the 
development of shear singularities in anisotropic scenarios (as has been shown in section IV). 

In the model (II. ip . regular sources of matter and radiation can remove the singularities 
thanks to the unconventional interplay between the matter and the geometry at very high 
energies. Due to the form of the gravity Lagrangian (II. ip . at low energies the theory recovers 
almost exactly the dynamics of GR because the connection coincides with the Levi-Civita 
connection of the metric up to completely negligible corrections of order ~ n'^p/Rp. At 
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Figure 5: Logarithmic representation of the shear as a function of K?p/Rp in radiation universes 
with different value of the anisotropy, which is controlled by the combination C^ = C^g + ^23 + ^' 
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In this representation, the difference between the curves is just a constant shift of magnitude log C 
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Figure 6: Evolution of the expansion as a function of K^p/Rp in radiation universes with low 
anisotropy, which is controlled by the combination C^ = CI2 + C'23 + ^31 ■ "^^^ case with C^ = 
corresponds to the isotropic flat case, 9"^ = 9H^. 



high energies, however, the departure is significant and that results in modified dynamics 
that resolves the singularity. The assumption that metric and connection are regarded 
as independent fields (Palatini variational principle) is at the root of this phenomenon, 
which could provide new insights on the properties of the quantum geometry and its 
interaction with matter. Because of this independence between metric and connection, the 
dynamics of our model turns out to be governed by second-order equations. As a result, the 
avoidance of the big bang singularity is not due to the existence of multiple new solutions 
of the field equations suitably tuned to get the desired result. Rather, the physically 
disconnected contracting and expanding solutions found in GR, which end or start in 
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Figure 7: Logarithmic representation of the expansion as a function of K?p/Rp in radiation uni- 



verses with high anisotropy, which is controlled by the combination C 
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singularities, are suitably deformed due to the non-linear dependence of the expansion on 
the matter/radiation density and produce a single regular branch (this non-linear density 
dependence is also manifest in loop quantum cosmology J7| and has recently been identified 
in [iTJ as a possible solution to the anisotropy problem). At low energies, the standard 
solutions of GR are smoothly recovered, and such solutions uniquely determine the high 
energy behavior. This should be contrasted with the same f{R,Q) Lagrangian formulated 
in the metric formalism, where due to the existence of additional degrees of freedom 
multiple new solutions arise and one must use an ad hoc procedure to single out those 
which recover a FRW expansion at late times. 

The fact that for all negative values of the parameter a one finds isotropic and anisotropic 
bouncing solutions in universes filled with dust indicates that with the Palatini modified 
dynamics the mere presence of matter is enough to significantly alter the geometry to avoid 
the singularity. Unlike in pure GR, there is no need for exotic sources of matter/energy 
with unusual interactions or unnatural equations of state. Regular matter is able by itself 
to generate repulsive gravity when a certain high energy scale is reached. And this occurs 
in a non-perturbative way. In fact, in the case of a radiation universe, for instance, the 
perturbative expansion of Q (see below equation (15.121) ) does not suggest the presence 
of any significant new effect as the scale Rp is approached. However, a glance at the 
exact expression (15. 9 p shows that there exists a maximum value for p, which is set by the 
positivity of the argument of the square root. Such limiting value is only apparent when 
the infinite series expansion of Q is explicitly considered. It is interesting to note that 
this type of non-perturbative effect arises in our theory without the need for introducing 
new dynamical degrees of freedom. In other approaches to non-singular cosmologies, the 
non-perturbative effects are introduced at the cost of adding an infinite number of derivative 
terms in the action (see 18| for a recent example). 

Additionally, since in radiation dominated universes {w = 1/3) the scalar curvature 
vanishes, R = 0, the mechanism responsible for the bounce in these models is directly 
connected with the Q = Ru^R^'^ term of the Lagrangian. In fact, all Lagrangians of the 



24 

form /(-R, Q) = f{R) + Q/Rp, will lead to the same cosmic dynamics^ ii w = 1/3, as is easy 
to see from the definitions fl2.6p - fl2.9l) and fl2.12p . This is a clear indication of the robustness 
of the models f{R, Q) = f{R) + Q/Rp against cosmic singularities. Note, in addition, 
that the anisotropic bounce always occurs when the maximum value of Q is reached, which 
emphasizes the crucial role of this term in the dynamics. 

On the other hand, the fact that this class of Palatini f{R, Q) actions can keep anisotropies 
under control for a very wide range of equations of state (including radiation and dust) 
without the need for introducing exotic sources (as in ekpyrotic models, which require 
w > 1), turns these theories into a particularly interesting alternative to non-singular 
inflationary models. 

To conclude, our investigation of anisotropies in Palatini f{R) and f{R,Q) models has 
been very fruitful. On the one hand we have been able to identify serious limitations of the 
f{R) models in anisotropic scenarios, namely, the existence of generic shear divergences, 
which makes these models unsuitable for the construction of fully viable alternatives to GR. 
On the other hand, we have shown that the model (II. ip is a good candidate to reach the 
goal of building a singularity free theory of gravity without adding new dynamical degrees 
of freedom. Whether this particular model can successfully remove singularities in more 
general spacetimes is a matter that will be studied in future works. 
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